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Abstract 

The propagation of low frequency seismo- acoustic waves in the Arctic Ocean ice 
canopy is examined through the analysis of hydrophone cind geophone data sets collected 
in 1987 at an ice Cctmp designated PRUDEX in the Beaufort Sea. 

Study of the geophone time series generated by under-ice explosive detonations reveals 
not only the expected longitudinad auid flexurad waves in the ice plate, but also an unex- 
pected horizontally-polarized transverse (SH) wave cirriving at a higher amplitude thcin 
the other wave types. The travel paths of aU three observed wave types are found to be 
refracted in the horizonted pleine cdong a line coincident with a known ridge sepairating 
the ice canopy locally into two distinct hedf- plates, the first of thin first yeax ice cind the 
second of thicker multi-yeeir ice. The origin of the SH wave appecirs to be near the detona- 
tion and not associated with the interaction of longitudined, flexural or waterborne waves 
with the ridge line. The need to determine the exact location of each detonation from the 
received time series highlights the dramatic superiority of geophones over hydrophones in 
this application, as does the ability to detect the cinomcdous SH waves and the refracted 
ray paths, neither of which axe visible in the hydrophone data. 

Inversion of the geophone data sets for the low frequency elastic p 2 LTcimeters of the ice 
is conducted initially by treating the ice as a single homogeneous isotropic plate to demon- 
strate the power of SAFARI numerical modeling in this application. A modified stationary 
phase approach is then used to extend SAFARI modeling to invert the data sets for the 
elastic pcirameters of the two ice half-plates simulteineously. The compressionad/shecir 
bulk wave speeds estimated in the half-plates, 3500/1750 m/s in the multi-yecir ice and 
3000/1590 m/s in the new ice, axe comp 2 Lrable to previously obtedned vedues; however, the 
compressionad/sheair attenuation vadues in the two hzdf-plates, 1.0/2.99 dB/A and 1.0/2.67 
dB/A, respectively, cire somewhat greater than previously measured values axid four times 
greater than estimates extrapolated from high frequency data. 

Thesis Supervisor: Dr. Henrik Schmidt 
Massachusetts Institute of Technology 
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Chapter 1 



Introduction 



The introduction describes the motivation and objectives for this thesis, and 
reviews its organization and content by chapter. 

1.1 Motivation 

The Arctic Ocean has grown increasingly important to national strategic interests 
in recent years; yet our understanding of the Arctic Ocean environment has lagged far 
behind that of the other major ocean systems. In particular, the modeling of low 
frequency acoustic propagation under the sea ice canopy in the Arctic Ocean has proven 
to be an elusive problem [1]. The difficulty has not been in general that the necessary 
tools to do this modeling are unavailable. For instance, Schmidt’s Fast Field algorithm, 
SAFARI [2], has proven to be a very capable package for solving propagation 
problems in a complex seismo-acoustic environment such as is presented by the deep 
Arctic. The difficulty has been that very little work has been done to obtain 
measurements of the starting parameters crucial to computing this propagation 
accurately, i.e., the elastic parameters of the ice canopy - compressional and shear wave 
bulk velocities and attenuation factors. As a result, previous modeling has been based 
largely on parameters measured in the laboratory or extrapolated from somewhat similar 
environments (freshwater lake ice, glacial ice, etc.). Recent investigations have 
suggested that parameters so obtained do not accurately reflect the Arctic environment. 
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and that measurements of the elastic parameters of the arctic ice are needed in situ 
before we can hope to model acoustic propagation in the Arctic Ocean accurately. 

1.2 Thesis Objectives 

The overall objective of this work is to improve our ability to accurately model 
acoustic propagation in the Arctic Ocean. To this end, data sets obtained in 1987 during 
seismo-acoustic propagation experiments conducted at an ice camp in the Beaufort Sea 
designated PRUDEX are studied extensively. The initial intent of this study was simply 
to apply advanced modeling techniques to the problem of inverting the hydrophone and 
geophone data as received at the PRUDEX arrays to obtain accurate measurements of 
the elastic parameters. Although inversion to obtain elastic parameters remains the 
focus of this work, unexpected phenomena observed in the propagation data have served 
to partially frustrate the immediate goal of obtaining the elastic parameters by increasing 
the difficulty of the inversion, while simultaneously contributing to the overall 
understanding of seismo-acoustic propagation in the Arctic by disclosing mechanisms 
in the propagation previously unobserved or unsuspected. In particular, arrivals 
characteristic of the refraction of all types of propagating waves at a linear discontinuity 
in the horizontal plane of the ice plate are presented. The author also attempts with 
apparent success to model this refraction and extend the inversion to obtain the elastic 
parameters of two separate kinds of ice cover, annual ice and multi-year ice, 
simultaneously. More importantly, completely unexpected horizontally polarized 
transverse (SH mode) waves are presented in the propagation data. Although the 
existing theories of seismo-acoustic propagation in an elastic plate have no mechanism 
by which an underwater explosion can generate SH waves in a sheet of ice floating over 
that explosion, these waves are present in the PRUDEX data sets at amplitudes greater 
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than any other wave type observed. These waves are important because their creation 
by out-of-plane scattering of other wave types may be a significant and heretofore 
unknown mechanism for the attenuation of acoustic energy entering the ice cover from 
the water. In addition to identifying their existence in the PRUDEX data, as a first step 
toward understanding their origin, the author investigates possible source locations for 
the SH waves relative the explosive source, the other wave types in the ice, and known 
discontinuities in the ice cover. 

1.3 Thesis Content 

Chapter 2 of this thesis lays the foundation for later work by reviewing the 
theory of wave propagation in a thin elastic plate under various conditions, focusing on 
the development of the three wave types commonly observed in such plates, the 
longitudinal plate wave, the flexural wave, and the transversely polarized SH wave. The 
second chapter then reviews briefly the numerical modeling tool used throughout this 
work, Schmidt’s SAFARI algorithm [3]. Chapter 2 concludes with a short discussion 
of the approach to the attenuation of elastic waves employed in this work. 

The third chapter introduces the reader to the experimental data as obtained at 
the PRUDEX ice camp. The background for the experiments is reviewed to make clear 
the need to determine through the analysis of seismo-acoustic propagation data some 
parameters which could have been measured precisely during the experiment. The 
nature of the acoustic source pulse used during the experiments, a key factor in later 
analysis, is described and explained using data received at a hydrophone array. Finally, 
the occurrence of the three principal wave types is identified in the data received at the 
experiment’s geophone array. 

Chapter 4 exists principally because the location and time of the underwater 
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explosive detonations used to excite seismo-acoustic propagation in the ice cover were 
not recorded relative to the receiving arrays during the PRUDEX experiments. In the 
process of determining shot time and location for use in later inversion, this chapter 
brings out the somewhat startling result that a simple array of four 3-axis geophones can 
be a much more effective tool for locating underwater sources than a larger hydrophone 
array in the water below the geophones. Geophone data is used in Chapter 4 not only 
to determine a much more accurate source location than is available using hydrophone 
data, but also to identify and analyze the refraction at the joint between two abutting ice 
half-plates of all wave types propagating in the horizontal plane of the ice sheet . 
Chapter 4 also reviews the evidence available to help identify the origin of the high 
amplitude SH waves which are visible only in the geophone data. 

Chapter 5 begins by explaining the fundamentals of the process of inverting 
response data for the elastic parameters of the propagating media, and then reviews 
previous work done to determine those parameters in arctic sea ice. The description and 
results of the inversion obtained by treating the ice canopy as a single homogeneous 
isotropic plate follow. These results serve to demonstrate the potential of SAFARI 
modeling of wave propagation; although based on the work of Chapter 4, the ice is more 
accurately modeled as two abutting half-plates with significantly different elastic 
parameters, and the results obtained by the single plate model are of questionable 
accuracy. To solve this problem. Chapter 5 introduces a method for using stationary 
phase analysis [4] in a somewhat modified form to extend two-dimensional SAFARI 
to model propagation of the flexural wave in the range dependent environment of the 
two abutting half-plates. The results of inversion using this modeling technique are then 
presented. 

Chapter 6 summarizes the results of Chapters 3, 4 and 5 and establishes their 
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importance relative to seismo-acoustic propagation in the Arctic Ocean. Chapter 6 also 
comments on additional work needed in this area. 
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Chapter 2 
Theory 

Prior to beginning the evaluation and analysis of any experimental data, it is 
important to have a firm underpinning in the theory behind that data. This chapter will 
review briefly the fundamental theory of propagation of elastic waves in an ice plate and 
the equations which characterize that propagation. It will then describe the principal 
tool used in this paper to solve these equations numerically, Schmidt’s SAFARI 
algorithms [3]. 

2.1 Propagation of Elastic Waves in a Plate 

An understanding of the unique characteristics of wave propagation in a thin elastic 
plate is essential to the study of seismo-acoustic waves in the arctic ice. Three 
fundamental wave types, longitudinal waves, flexural waves and horizontally-polarized 
transverse (SH) waves are commonly observed propagating in floating ice sheets 
[5][6]. It is useful to look first at the origins of the three wave types in a free elastic 
plate and then extend those results to a plate bounded by a liquid half-space on one side 
and a vacuum (or air) on the other. 

2.1.1 Elastic Waves in a Free Plate 

Consider waves propagating in the positive x direction in a laterally infinite 
homogeneous isotropic plate bounded on both sides by vacuum as shown in Figure 2-1. 
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Figure 2-1: Two dimensional geometry for infinite plate of thickness 2h 
bounded by vacuum above and vacuum or liquid below. 

If displacements in the plate satisfy the Navier equation, 

pM = /+(A.+2p)V(V«)- |iVx(Vx«) , 



where X and p are the Lame constants, then Lame’s theorem states [4] that a scalar 
potential <j) and a vector potential exist which satisfy 



u = V4> + Vx\l; 


(2.2) 


V*i|r = 0 


(2.3) 




(2.4) 


-e-i 

II 

^.1- 


(2.5) 
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where 



^2 ^ X+2\i 
P 



( 2 . 6 ) 



P 



2 



P 

P 



(2.7) 



As a trial solution, consider plane wave solutions of the form 

$ = 

Substitution of these solutions into (2.2) through (2.5) yields the following set of four 
equations in eight unknowns for the potentials, 

4) = (2.10) 

\|r = (2.11) 

Tjf = (2.12) 



tjr = , (2.13) 

The wave nature of the potentials, and thus the motion, is clearly seen in these 
equations. 

(2.10) through (2.13) can be used to develop seismo-acoustic propagation in any 
layered (two dimensional) environment, and various types of body and surface waves 
arise depending upon the media and the boundary conditions. The quantities a and P 
are the compressional (P-wave) and shear (SV/SH-waves) bulk velocities; Iq is the 
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horizontal wavenumber and and are the vertical wavenumbers for compressional 
and shear waves. No attempt is made in this paper to review the theory supporting any 
propagation other than that occurring in a thin elastic plate; although references to other 
wave types are used to relate propagation in the thin elastic plate to other results for 
readers familiar with seismo-acoustics. 

From this point, one of the most complete developments of the basic 
characteristics of waves in a free plate is provided by Graff based on work by Meeker 
and Meitzler [7]. For simplicity in following this development, (2.10) through (2.13) 
are rewritten in terms of sines and cosines (with different constants). 



<f) = (i4cosA:^^z + (2.14) 

ilfj^ = (CcosA^pZ + DsinA:^pZ)e^***'“*^ (2.15) 

tjfj, = (£dosfe^pZ + Fsm^^pz)e^*''"*^^ (2.16) 

t|/j = (Gcos^^pZ + /fsinA:^pZ)c*^*'*"“‘^ . (2.17) 



If (2.2) is simplified for dependence only on x and z in the two dimensional problem 
under consideration, and the potentials of (2.14) through (2.17) are substituted into the 
result, then particle displacements are given by 

6(J) 3i|;^ 

(2-18) 

= [ikJiAc(xK^^z+Bsmk^^z)+k^pi-E^iiik^pZ+Fcosk^pZ)]e‘^''‘^''^‘^ 
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y 



(2.19) 



dz 3x 

= [ (- CsiiiJr pZ +Dcos^^pz) +iA:/Gcos^^pZ +/TsinA:^pz)]e 



a* 9i|iy 

U = -2- L 

" 3z Bx 

= [^^^(-i4sinisr^^z+-Bcos^^^z)-i^,(£tosfc^pZ+Fsin^^pZ)]c‘^*** . 



( 2 . 20 ) 



The basic nature of wave motion in elastic solids emerges from (2.18) through (2.20) 
in the decoupling of transverse (SH) motion from radial and vertical motion (SV/P); Uy 
depends on X|/^ and X|/^, while u^ and u^ depend on \|/y and <{>. 

The generalized form of Hooke’s law for the stress tensor, Xy, in a linearly 
elastic solid, 






( 2 . 21 ) 



where the strain tensor, ey, is given by 



1 






( 2 . 22 ) 



reduces in the homogeneous isotropic plate to 



du^ 0a, 

T_ = (A.+2p ) — - + X 

= 0z dx 



= P 



^ du 0a ^ 
— - + — - 
^ dx dz 
du„ 



(2.23) 



X = u — ^ . 

^ ^ 0(y 

On the upper and lower faces of the plate (z=±h) the boundary conditions 



between plate and vacuum are simply 



(2.24) 






Equations (2.18) through (2.20), (2.23) and (2.24) combine to form 6 equations 
in the eight unknowns A, B, C, D, E, F, G and H. Recalling that Lame’s theorem also 
provides that the divergence of the vector potential is zero (2.3), in the simplified 
geometry of the plate this relation becomes 



dx dz 



= 0 



(2.25) 



Substituting (2.15) and (2.17) into (2.25) yields two additional equations, for a total of 
eight equations in the eight unknowns. In matrix form this system of 8 equations is 
written as 
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G 
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= 0 



(2.26) 



where 



- 20 - 



b = 

c = 4 

J = x^Ap 

e 

/ = ^.'-4 




c„ = coskh 

OL ZCL 

s„ = sink,„h 

OL ZCL 



Cp 



cosk^ph 



Sp = sinA^p/t . 



(2.27) 



Using a number of straightforward matrix manipulations, it is easy to show that the 
system of (2.26) is equivalent to 
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(2.28) 



A necessary and sufficient condition for the existence of a solution to the homogeneous 
system Ax=0 is that det(A)=0 [8]. The format of (2.28) is advantageous because the 
determinant of that special form of an 8x8 matrix expands to the simple form 



as^ -bsp 




-ds^ csp 


/Cp 




-hs^ gs^ 



/Sp 




CCp dcp 


ac^ fccp 




gCp /xcp 



0 . (2.29) 



where coefficients B and E are associated only with the elements of the first, G and D 
the second, A and F the third, and C and H the fourth determinant in (2.29). Thus four 
separate solutions to (2.28) must exist, such that only the pair of coefficients associated 
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with each 2x2 determinant is non-zero, while the other six coefficients vanish. 

Looking first at the coefficient pairs A/F and B/E, if only A,F?K), then 
substitution in (2.18) through (2.20) yields 

My = 0 (2.30) 

and if only B,E?K), the same substitution produces 

My = 0 (2.31) 

Ux = (kx^Bcoskx^z-ikxEcoskxpz)e‘^^^ . 

Solutions (2.30) and (2.31) are polarized in the longitudinal and vertical directions only; 
(2.30) is symmetric with respect to the x/y plane while (2.31) is antisymmetric with 
respect to the x/y plane. Additionally, setting the 2x2 determinant in (2.29) associated 
with each solution’s coefficient pair to zero generates an equation which describes the 
required relation between vertical and horizontal wave numbers for that solution to exist, 
i.e., the frequency equation for that solution. Setting the 2x2 determinant for 
coefficients A and F to zero yields 

-es^bc^ -fs^ac^ = 0 , (2.32) 

which can be rewritten as 



tan^^p^ ^ 

and since 



(2.33) 
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where 






1 



}:^Y_ 

1 






1 



-kl-lk] 

p2 

1 



,2 ,2 

K^-K 



(2.34) 



kl = k; + kl^ and kl = k] + kl^ , 



(2.35) 



(2.33) becomes 



~ ' (k-k^^f 



(2.36) 



(2.36) is the well known Rayleigh-Lamb frequency equation for symmetric plate waves. 
In the antisymmetric case (coefficients B,E?tO) the determinant yields similarly 



tan^zp^ ^ __(V^P^ 
tan^r-^ Ak^k k 

za ‘t/tjp /t^p 



(2.37) 



Although (2.36) and (2.37) were originally obtained slightly more than 100 years 
ago by Rayleigh and Lamb, the associated frequency spectrum was not completely 
understood until much later. Figure 2-2 illustrates the complete free plate frequency 
spectrum as obtained by Mindlin [9] for a value of Poisson’s ratio of d= 0.31. In 
Figure 2-2 the (positive or negative) real and imaginary values of k,^=(2h/jc)k,j are 
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imaginary _ real 

k^=(2h/7i)k. 



Figure 2-2: Frequency spectrum, kj=(2h/7:)kj versus kp=(2h/7t)kp, for \)=0.31, 
showing symmetric (thick lines) and anti-symmetric(dotted lines) modes (from 
Mindlin [9]). 

plotted against kp=(2h/7i)kp. While the periodic nature of the tangent function 
introduces an infinity of modes. Figure 2-2 also shows that except for the lowest order 
fundamental symmetric and anti-symmetric modes, all modes have a cutoff frequency 
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below which the propagating wavenumbers are imaginary and waves are evanescent in 
the horizontal direction. The cutoff frequencies for the first symmetric and 
antisymmetric modes above the fundamental ones are given by 




and 




(2.38) 



respectively [7]. Treating a typical arctic ice sheet (h=3m, a=3500m/s, P=18(X)m/s [14]) 
as a free plate, below about 300 Hz only the fundamental antisymmetric and symmetric 
modes can propagate. Further, for the case in which the wavelength is long compared 
to the thickness of the plate (thin plate assumption), i.e., k<jh-40 and kph^O, (2.36) can 
be reduced easily to a simple expression for the fundamental symmetric mode by 
replacing the tangent functions with their arguments 



^ =4 
p2 



1-Pi 



(2.39) 



where c,p is the phase speed for the wave. The lowest order symmetric mode in a thin 
plate is called the longitudinal or plate wave. The longitudinal wave is non-dispersive 
in the thin plate limit since (2.39) indicates that the phase speed is not dependent on 
frequency. A similar expression can be obtained for the lowest order mode of the 
antisymmetric case in the thin plate limit. By replacing the tangent function with its 
Taylor series and discarding higher order terms, (2.37) reduces to 






2 \ 



i-i- 



a" 



(2.40) 



3 

Here Cf is the phase speed for the lowest order antisymmetric wave, generally referred 
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Figure 2-3: A symmetric (longitudinal) wave in a free plate seen as the 
superposition of a pair of P waves and a pair of SV waves incident on the faces 
of the plate (after Redwood [1 1]). 

to as the flexural wave. The flexural wave is clearly a dispersive wave. As implied by 
the form of (2.30) and (2.31), longitudinal and flexural waves in a free plate may be 
thought of as constructive interference of multiply-reflected P and SV body waves 
[10]. Figure 2-3 (after Redwood [11]) demonstrates this point of view for a 
longitudinal wave. 

Looking now at the remaining two coefficient pairs in (2.28), if only C,H?tO, then 
substitution in (2.18) through (2.20) produces 

< 2 . 41 ) 

and if only D,GvO, then the same substitution produces 
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(2.42) 



X z 

Uy = [i-k^pD + ik^G)cosk^pZ]e^‘'^"^'^ 

These two solutions are transversely polarized SH waves; (2.41) is antisymmetric with 
respect to the x/z plane, and (2.42) symmetric with respect to the x/z plane. The 
frequency equation from the 2x2 determinant in (2.29) associated with the antisymmetric 
case is 

cCpAcp-JcpgCp = 0 , (2.43) 

which simplifies to 

Solutions to (2.44) exist only for 

k^ph = , n=0,l^,3... . (2.45) 

2 

Similarly, for the symmetric case 

= (J » (2.46) 

for which solutions exist only for 

= m It, ffi=0, 1,2,3,... . (2.47) 

If the results of (2.45) and (2.47) are combined, then SH mode solutions exist for 

Jt./i = — , n=0,l,2,3... , (2.48) 

zp 2 

with n odd the antisymmetric and n even the antisymmetric modes. (2.48) can be 
shifted to a more revealing form by writing it in terms of kp and k„ 
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Figure 2-4: Frequency spectrum, versus kp, for the symmetric(-) and 
antisymmetric(— ) SH modes of an infinite free plate (after Graff [7]). 



ik^h) = 



n% 



\2 



<khf 



(2.49) 



(2.49) is plotted in Figure 2-4 on axes again scaled by 2h/7i. Note in Figure 2-4 that all 
but the fundamental symmetric mode and all antisymmetric modes have cutoff 
frequencies given by 



o> 




(2.50) 



below which modes become evanescent (imaginary horizontal wavenumber); however, 
the fundamental symmetric mode, n=0 in (2.49), is a nondispersive propagating wave 
independent of frequency. Since the cutoff frequency given by (2.50) is relatively high 
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for a typical arctic ice sheet (again, ~300 Hz in 3 meter ice), it is likely that only the 
fundamental symmetric SH mode will be of concern in this investigation. 



2.1.2 Elastic Waves in a Floating Plate 

The problem of a thin elastic plate over a liquid half-space differs from the free 
plate problem only in the boundary conditions at z=h. In place of (2.31), the boundary 
conditions are unchanged at the upper boundary with vacuum. 



«i 



= T 



«i 



= t 






= 0, z=-h 



but at the lower boundary they become 



(2.51) 
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= 0 . 
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«. = W 






z=+/i 



(2.52) 



as is appropriate at a liquid/solid interface [7]. The subscript 1 refers to the plate and 
the subscript 2 the underlying liquid half space in (2.51), (2.52), and the work which 
follows. 

Unfortunately, when the floating plate boundary conditions are used with 
equations (2.18) through (2.20) and (2.23), the new system of nine equations in nine 
unknowns cannot be analyzed as readily as the free plate system. The P and SV wave 
motions in the floating plate no longer reduce to purely symmetric and antisymmetric 
modes due to the presence of the liquid. Press and Ewing [12] studied the P and SV 
modes by making the simplifying assumption that A,=p, and derived an exact expression 
for the period equation. 



where 



P(2Q+5coshVjft coshvi/i)+Q6siiihVj/i sinhvj/i = 0 , 



(2.53) 
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5 = pjCtjCvf 



PiP^a 

V = , v' = ik^p , (2.54) 

P = (vj^+fc^^coshvj/i sinhvjA - 4vjVjj!:^siiihVj/i coshvj/i , 

Q = (vj^+t^^sinhVjA coshvjA - 4v jVjit^coshv^ft sinhv j/i . 

Press and Ewing were further able to show that for long waves in a thin plate (k^h 
small), flexural waves analogous to the free plate’s antisymmetric case exist and 
propagate with a period equation given by the approximation 




1 -— 
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3pj ^ ^ 
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(2.55) 



and analogous longitudinal (symmetric case) waves exist and propagate with a period 
equation given by the approximation 



= 
Pi 
/ 

1- 



PI 



2\ 1 



a, 






3P1P1 



-2Nir 



®2P2 



,-ii 

2 



a 



/ L 



1 - 



4P? 



a, 



i-ii 

2 






a 



(2.56) 



Note in particular that the real part of the longitudinal wave velocity is unchanged from 
the free plate (2.39), but that the liquid half-space adds an imaginary part which 
represents attenuation proportional to the wavenumber cubed. In the short wavelength 
limit (k^h large). Press and Ewing showed that (2.53) produces Rayleigh waves on the 
free surface and Rayleigh and Stoneley (Scholte) waves at the ice/water interface. These 
interface waves are not important for the ice thicknesses and frequencies of concern in 
this work and will not be discussed further. 
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The change in boundary conditions is irrelevant to the SH wave, the derivation 
of which involves neither nor u^. The SH wave propagates unchanged in a plate 
regardless of the fluid (or vacuum) on the faces of the plate. Thus the solution for the 
SH modes in the floating plate is exactly that for the free plate discussed in the previous 
section. 

2.2 Numerical Solutions 

Numerical solutions to the transcendental characteristic equations for the floating 
plate have been obtained generally by making simplifying assumptions. As seen in the 
previous section, Press and Ewing assumed that the Lame constants were equal and 
looked at the solution in the limit as the wavelength became very small or very large. 
A number of thin plate theories have been studied which model a fluid loaded plate in 
such a way as to account for only the lowest flexural and possibly longitudinal modes. 
Although Langley [13] has shown that the thin plate approximations can provide 
accurate results below the cutoff frequencies for the higher order symmetric and 
antisymmetric modes, with modem computing facilities and available tools these 
approximations are no longer necessary to achieve quick and accurate results from the 
exact equations. Stein [14] used a specialized computational routine to solve the P- 
SV system of equations for the floating plate numerically; however Schmidt [2] has 
developed a much more flexible tool to apply to this problem, the Seismo- Acoustic Fast 
Field Algorithm for Range-Independent environments (SAFARI). 

2.2.1 Full Wavefield Global Matrix Solution 

The SAFARI approach to solution of seismo-acoustic propagation in a 
horizontally-layered environment is at its heart an expansion of the techniques of solving 
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the depth-separated wave equation originally applied to acoustic propagation by Pekeris 
[15] and later extended to seismic propagation by Ewing, Jardetzky and Press [10]. 
To see how this technique works in SAFARI [3], the wave equation for compressional 
waves (2.5) is rewritten for a single horizontal layer in cylindrical coordinates with the 
assumption that the environment is axi-symmetric, range-independent (i.e., two 
dimensional) and all sources are on the z-axis, 









W^,t) =/,(z,f)5(r) 



If the Fourier transform. 



(2.57) 



F(o>) * , 



(2.58) 



and the Hankel transform. 



G(k) = f'g(rWo(k/)rdr , 

JQ 



(2.59) 



are both applied to (2.57), the result is an ordinary differential equation in z only. 






[dz^ 






(2.60) 



/ 



2n 



where Iq is the horizontal wavenumber and k„=a)/a as before. The solution to (2.60) 
is the depth-dependent Green’s function, given at some radial frequency co by 



= 



(2.61) 



where Op is some particular solution to (2.60), O' and O^ are two independent solutions 
to the homogeneous form of (2.60), and A' and A^ are coefficients determined by the 
boundary conditions. If the z-dependence of Iq, is restricted to cases for which (2.60) 
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can be solved analytically, then solution of (2.60) in a multi-layer system amounts to 
solving for the arbitrary coefficients of (2.61). If a is taken as constant and there are 
no sources, for instance, then solving for the depth-dependent Green’s function yields 

+ A*e^ , (2.62) 

and the potential is 



<|)(r,z) = j^A~e~*^ + A*e^\lQ(kr')k^^ . (2.63) 

For solid media, SAFARI takes P as constant (in each layer) and solves for 
^F(r,z) in a manner similar to that above; although as implemented in SAFARI the 
vector \|/(r,z) is replaced by (in cylindrical coordinates) x/ (r,z)=3\|/e/3r. The equations 
for the displacements and stresses become 



and 



az r dr dr 



(2.64) 
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3z dr 



Thus the depth-dependent Green’s function is 






and the appropriate solution for the \j/ potential is 



(2.65) 



( 2 . 66 ) 
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(2.67) 



iJfVrZ) = + 



Substituting (2.63) and (2.67) into (2.64) and (2.65) produces a set of four 
equations in the four unknown coefficients for the displacements. 






( 2 . 68 ) 



and the stresses. 



Tjr,z)= iif^(2k^,-ki)(A-e-^ 

-lik^^^iB-e'^^^^ - B^e‘^^^]f^(k/)k/ik^ 
x^(r,z) = \ij^ikJc^^(A-e'^ -A*e^) 



(2.69) 



This set of equations must be satisfied on both sides of each boundary between layers 
such that the boundary conditions appropriate to the adjacent layers 

(liquid/solid/vacuum) are satisfied. Since the boundary conditions must be satisfied at 
all ranges, the kernels of the integrals in (2.68) and (2.69) must also satisfy the 
boundary conditions, leading immediately to a linear system of equations in the 
unknown coefficients. A', A'^, B- and B"^, at every horizontal wave number k,. 

In SAFARI the linear system of coefficients is solved numerically to determine 
the depth-dependent Green’s functions at desired depths for a discrete set of 
wavenumbers; the set of wavenumbers chosen must be sufficient to allow the numerical 
determination of the inverse Hankel transforms in (2.68) and (2.69) for the desired 
ranges from source to receiver. The resulting transfer functions are determined at a 
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discrete set of frequencies, and again the set of frequencies chosen must be sufficient 
to support the numerical calculation of the inverse Fourier transform to determine the 
total response in the time domain. In principle, the SAFARI solution is exact, limited 
only by the accuracy of the numerical methods used to solve the linear systems and 
determine the inverse wavenumber and frequency transforms, and the need to describe 
the environment by discrete layers within which the wave equation is separable. 
Clearly, SAFARI must be applied with care and knowledge, in particular weighing the 
requirements of wavenumber and frequency sampling against the limitations of machine 
memory and processing time. For an investigator with a firm background in the 
fundamentals of wave propagation and numerical analysis, and with a little experience 
with the code, SAFARI is an ideal tool to apply to the problem of wave propagation in 
an elastic plate. 

2.2.2 Attenuation 

The treatment of elastic waves in a plate to this point has assumed frictionless 
propagation, but in reality some energy is lost due to internal friction with each cycle 
of stress. As developed in Aki and Richards [4], a dimensionless measure of this 
friction is 



1 



a£ 



(2.70) 



Q(u) 2nE 

where E is the peak strain energy and -aE is the energy lost in each cycle. If 
attenuation is assumed to be a linear phenomenon, wave amplitude. A, is proportional 
to in a linearly elastic solid. If Q»1 is also assumed, (2.70) can be rewritten 



1 ^ _J_±i 

0(0)) n A 



(2.71) 
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For attenuation of a wave propagating in the x direction 




and (2.71) becomes 



(2.72) 



^ = __5d_ = 
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for which the solution is 



(2.73) 



A(x) = AqC 




(2.74) 



This relation can be viewed as treating the wavenumber of a propagating wave as a 
complex value, i.e.. 




(2.75) 



The explicit assumptions of linearity used in the above development are made 
in the SAFARI code, and this approach to attenuation is adopted in this paper. It is 
common in ocean acoustics to express attenuation in dBA, so that linear frequency- 
dependent attenuation is given by the parameter y, where 



^ 207ilogg (2.76) 

Q 

This convention is used in SAFARI and throughout this paper. 
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Chapter 3 

Experimental Measurements 

Before embarking upon detailed discussions of source location and inversion in 
later chapters, it is useful first to review the nature of the experimental measurements, 
both to understand the motivation for the chosen approach to propagation analysis 
(including the need for extensive work with localization of experimental sources) and 
to observe certain aspects of the measured wave forms which are in themselves 
astonishing. This chapter will first review the seismo-acoustic propagation experiments 
which generated the data used in later chapters and then will take a close look at the 
data collected, including the time series of longitudinal, flexural and SH mode waves 
generated by under-ice, underwater explosive charges. 

3.1 The Experiments 

All data sets included in this investigation were collected during experiments 
conducted during March and April 1987 at an ice camp designated PRUDEX, located 
in the Beaufort Sea approximately 100 nautical miles north of Prudhoe Bay, Alaska. 
The PRUDEX ice camp was established as a joint effort of the Woods Hole 
Oceanographic Institution, the Massachusetts Institute of Technology and the Polar 
Science Center of the University of Washington Applied Physics Laboratory with the 
primary objective of testing an Arctic Remote Autonomous Measurement Platform 
(ARAMP) data buoy [16]; however, other experiments were included in the program. 
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Figure 3-1: PRUDEX hydrophone and geophone array layout on the x/y plane 
used throughout this paper, hydrophones H0-H15 suspended at a 60 meter water 
depth, geophones G1-G4 frozen into the ice. 

Originally intended only as engineering experiments designed to verify techniques and 
procedures for full scale experiments in later years, seismo-acoustic ice propagation 
experiments were conducted to observe the response of the PRUDEX ice canopy to 
underwater detonations. Preliminary review of the propagation data sets recorded during 
PRUDEX indicated that the data sets were of such high quality as to warrant a much 
more detailed investigation than had been planned when the sets were obtained. 

Two co-located receiving arrays were used to record seismo-acoustic data at the 
PRUDEX camp. As shown in Figure 3-1, waterborne waves were received at a sixteen 
element hydrophone array. Each hydrophone was suspended at a depth of 60 meters at 
the locations shown. No exact hydrophone position data were available during the data 
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recordings of interest. Hydrophone response was sampled at 1000 Hz and recorded in 
digital form on an optical disk [17]. Waves propagating in the ice were received by 
an array of four 3-axis geophones. Geophones were frozen in the ice at the locations 
shown in Figure 3-1. The geophone data were recorded on wideband magnetic tape and 
later sampled at 1000 Hz and digitized by the author using the Woods Hole 
Oceanographic Institution’s MIZEX Analog to Digital Converter [18]. 

Seismo-acoustic propagation data sets were generated during PRUDEX with a 
series of eight sets of underwater explosive detonations (shots), each set conducted in 
the early morning (Universal Time) on successive days in late March and early April. 
Each set of shots was conducted at a different location at ranges of from 300 to nearly 
600 meters from the receiving arrays. Each set consisted generally of six separate shots 
conducted at source depths of 2, 4, 6, 8, 16 and 32 feet below the ice surface using 
various amounts of explosive charge. Table 3-1 summarizes the experimental shots 
conducted at the PRUDEX ice camp for which data were available. Of all shots 
completed, two shots in each of four data sets proved useful for the propagation analysis 
of this paper. Useful shots were primarily those which excited strong flexural and 
longitudinal waves in the ice, as well as producing a measurable geophone response to 
the waterborne acoustic wave generated by the detonation as the wave passed the ice 
beneath the geophone. The flexural waves were essential to the propagation analysis 
and inversion, while the longitudinal waves and the waterborne acoustic wave responses 
proved to be necessary for determining shot location. All of the required responses were 
generated by detonations at depths of 8 and 16 feet, at ranges of 500-600 meters, and 
using 1 or 2 feet of detonating cord (primacord) as the explosive source. 

Because the propagation data sets were planned only to verify equipment and 
procedures and not with subsequent analysis in mind, much of the supporting 
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date 


shot 

series 


shot desig. 


time (U.T.) 


depth (feet) 


charge 


3/26/87 


B 


B2 


0537:41 




6 drams 


B3 


0542:03 




1 ft cord 


B4 


0546:14 




2 ft cord 


3/27/87 


c 


Cl 


0418:37 


♦ 


4.5 drams 


C2 


0419:58 




6 drams 


C3 


0423:29 


♦ 


1 ft cord 


C4 


0427:24 




2 ft cord 


C5 


0436:00 


★ 


4.5 drams 


C6 


0438:51 




6 drams 


C7 


0442:36 




1 ft cord 


C8 


0445:33 




2 ft cord 


C9 


0456:34 


♦ 


4.5 drams 


CIO 


0458:06 


♦ 


6 drams 


Cll 


0501:54 




1 ft cord 


C12 


0504:29 




2 ft cord 


3/29/87 


D 


D1 


0534:56 




4.5 drams 


D2 


0537:40 


♦ 


6 drams 


D3 


0543:07 




1 ft cord 


D4 


0547:14 




2 ft cord 


3/31/87 


F 


F3 


0546:21 




1 ft cord 


F4 


0550:48 


♦ 


2 ft cord 


4/1/87 


G 


G1 


0442:06 


2 


2 ft cord 


G2 


0445:26 


4 


2 ft cord 


G3 


0450:45 


8 


2 ft cord 


G4 


0454:46 


16 


2 ft cord 


4/2/87 


H 


HI 


0548:25 


64 


2 ft cord 


H2 


0553:28 


32 


2 ft cord 


H3 


0557:14 


16 


2 ft cord 


H4 


0601:15 


8 


2 ft cord 


H5 


0604:39 


4 


2 ft cord 


H6 


0606:44 


2 


2 ft cord 


4/3/87 


I 


11 


0558:05 


64 


2 ft cord 


12 


0602:32 


32 


2 ft cord 


13 


0607:33 


16 


2 ft cord 


14 


0611:42 


8 


2 ft cord 


15 


0614:32 


4 


2 ft cord 


16 


0617:12 


2 


2 ft cord 



Table 3-1: Summary of experimental shots recorded during PRUDEX propagation 
experiments. (*) indicates no depth recorded on shot log, "cord" refers to primacord, 
"drams" to explosive weight. 
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information which could have greatly simplified analysis was not recorded. Principal 
among the deficiencies in recorded data was that, except in only the most general terms, 
neither shot location nor shot time was recorded relative to the receiving arrays. As a 
result, a mandatory and non-trivial preliminary to detailed propagation analysis was 
locating each shot using only the seismo-acoustic data received at the arrays. A second 
deficiency was the failure to conduct a geophone calibration either before, during or 
after the experiment. Additionally, because the data from the two different arrays was 
recorded on completely different systems and no attempt was made to align the two 
systems precisely, the hydrophone time series could not be synchronized with the 
geophone time series before analysis began. 

Further complicating the propagation analysis were known discontinuities in the 
ice ctmopy at PRUDEX; the explosive shots were made beneath the ice camp’s runway, 
a relatively thin floe of new ice, while the receiving arrays were located on or below an 
abutting floe of thicker multi-year ice. Again, as detailed analysis was not anticipated, 
the surface geometry of the floe abutment was not surveyed, nor were any useful ice 
thickness or under-ice surveys conducted. In fact, the only available ice thickness data 
were recollections of the personnel who drilled ice holes in support of hydrophone array 
and under-ice shot placement. The runway was informally reported as being about one 
meter thick, and the array ice reported as being somewhat variable, about three meters 
in thickness. 



3.2 The Observations 

The purpose of this section is to review the time series observed during the 
PRUDEX propagation experiments and associate observed wave forms with the 
appropriate wave types discussed in Chapter 2. Additional characteristics of the 
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rm 





time after detonation (sec) 

Figure 3-2: Waterborne acoustic waves as recorded at the output of PRUDEX 
array hydrophones (from top to bottom) #8, #3 and #0 in response to 
experimental under-ice shot number F3. 

observed time series which greatly simplify analysis and inversion are identified and 
explained. 

3.2.1 Hydrophone Data 

Figure 3-2 shows a typical time series resulting from the detonation of 1 foot of 
primacord under the ice at a water depth of 8 feet and a range from array center of 570 
meters, as recorded at the output of several PRUDEX array hydrophones. A striking 
feature of each of the traces in Figure 3-2 is the characteristic shock wave and bubble 
pulse pressure signature of an underwater explosive detonation. Bubble pulses, the 
pressure peaks which follow the initial detonation shock wave, result from successive 
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oscillations of the globular mass of gaseous materials that remains after the detonation 
is completed, each successive pulse being weaker than the preceding one as remaining 
energy is dissipated. Generally, only the first several bubble pulses are strong enough 
to be observable. Although the peak pressure of the first bubble pulse is about 40% that 
of the shock wave, at lower frequencies the energy density is actually higher in the first 
bubble pulse than in the shock wave. The principal peak in the combined spectrum of 
an underwater detonation occurs at a frequency of 1/T, with T being the interval 
between shock wave and first bubble pulse [19]. As can be seen in Figure 3-2, after 
low pass filtering to prevent aliasing in the data acquisition system, the first bubble 
pulse is actually higher in amplitude than the shock wave due to its higher energy 
density at low frequencies. This characteristic is common to most shots analyzed. 

3.2.2 Geophone Data 

Geophones in the PRUDEX array were implanted in the ice with principal axes 




Figure 3-3: PRUDEX geophone array layout showing alignment of principal 
axes on each geophone. 
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aligned as shown in Figure 3-3; however, a much more useful aspect is obtained by 
resolving the output of the x and y axis geophones to axes corresponding to the radial 
and transverse directions relative to the propagation path of the appropriate shot. 
Figure 3-4 shows the time series corresponding to the same shot shown in Figure 3-2 
as recorded at one 3-axis geophone and resolved to the direction of propagation. 

The first arrival on the radial particle velocity trace in Figure 3-4 is the 
longitudinal plate wave, the pulse corresponding to the detonation shock wave arriving 
first followed by a larger one corresponding to the first bubble pulse. As the 
longitudinal wave is only slightly dispersive, the characteristics of the underwater 
detonation are retained in its time series; unfortunately later multiple arrivals generated 
due to complex interactions with discontinuities in the ice tend to confuse the pattern 
for the subsequent weaker bubble pulses. The longitudinal wave is also visible in the 
vertical geophone, although much less so than in the radial geophone, and later multiple 
arrivals dominate the time series until the flexural wave begins. 

The flexural wave is the strongly dispersive wave beginning about 0.7 seconds 
after detonation in Figure 3-4 and continuing to the end of the trace. The flexural wave 
dominates the vertical geophone output, and is clearly visible in the radial geophone 
output. The dispersive nature of the flexural wave destroys the shock wave/bubble pulse 
characteristic of the response; however comparison with shots made at the same location 
but at different depths, hence with different bubble pulse intervals, shows that the 
flexural wave behaves as though its time of origin is the time of the first bubble pulse, 
not the initial shock. This correction is significant when determining the phase and 
group velocities of the flexural wave for comparison with theoretical results. 

A second distinct set of pulses on the vertical geophone, arriving at about 0.4 
seconds after detonation in Figure 3-4, correspond to the response of the ice plate by the 
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Figure 3-4: Time series as observed in the (top to bottom) radial, transverse and 
vertical directions on PRUDEX geophone #4 in response to experimental shot F3. 



- 45 - 



time after detonation (sec) 



waterborne acoustic pulse arriving at the underside of the floe directly below the 
geophones. Since the travel path in the ice is short, the pulses replicate nearly exactly 
the characteristics of the hydrophone arrivals in Figure 3-2. At the liquid/solid interface 
the energy transmission is only in the direction of the normal to the surface, and the 
waterborne pulse is seen principally in the vertical geophone; although discontinuities 
in the underside of the ice produce a very small response in radial and transverse 
geophones as well. The responses to the waterborne pulse carry no information useful 
directly in the inversion for the elastic parameters of the ice; nonetheless, it will be 
shown in Chapter 4 that these pulses are indispensable to accurately determining the 
location of the experimental shots relative to the array. 

Figure 3-5 contains two expanded plots for the vertical axis output of one 
geophone, with the response from two separate shots made at the same location overlaid 
for comparison. One shot was made using 1 foot of detonating cord at a depth of eight 
feet, and the other shot using 2 feet of cord at the same depth. The first plot in 
Figure 3-5 shows in detail the geophone response to the waterborne wave as it passes 
beneath the ice on which the geophone rests, as well as the change in the bubble pulse 
interval due to the change in the size of the charge. The second plot is of the flexural 
wave and demonstrates both the repeatability of the wave and the slight offset in the 
onset of the flexural wave introduced by the difference in the bubble pulse intervals. 

Contrary to all expectations based on the theory of plate wave propagation, the 
largest amplitude response to an underwater detonation in each of the 3-axis geophones 
occurs only on the transverse geophone at a velocity corresponding roughly to the 
expected shear velocity in the ice. Further, the wave form of this transverse response 
retains generally the shock wave/bubble pulse characteristic of the source, indicating that 
this arrival is at most only slightly dispersive. The inescapable conclusion to be drawn 
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Figure 3-5: Vertical geophone #4 response for experimental shots F3(-*) and 
F4(-); top, the response to the waterborne acoustic wave as it passes under the 
ice, and bottom , the flexural wave. 



from this combination of characteristics is that the largest amplitude wave propagating 
in the ice as a result of an under-ice, underwater detonation is the fundamental 
symmetric SH mode! This result is surprising because, as discussed in Chapter 2, there 
should be no interaction between SH waves in the plate and acoustic waves in the water, 
thus, an underwater detonation should not excite SH waves. Additionally, in a 
homogeneous, isotropic plate there is no mechanism to convert longitudinal and flexural 
(SV and P) waves into SH waves. Either the detonation does excite SH waves in a 
manner not understood, or out-of-plane scattering at discontinuities in the ice sheet 
provides a significant coupling between SH and SV/P modes. Either conclusion would 
be significant, in that most work in the arctic environment to date has ignored the SH 
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Figure 3-6: Hodograph for the x and y axes of geophone #3 showing the 
response to shot F4 at times from 0.1 to 0.323 seconds(-) and 0.323 to 
0.5 seconds(-*) after detonation. 

mode in the ice as insignificant relative to ocean acoustics. The assumption that the SH 
mode may be safely discarded appears to be questionable based upon these observations. 

Figure 3-6 shows a typical hodograph constructed from x and y geophone outputs 
during an under-ice shot. The hodograph further illustrates the clear radial polarization 
of the longitudinal wave arrival, followed by the equally clear transverse polarization 
of the later SH wave arrival. Similar geophone hodographs have been constructed 
during preliminary analysis of under-ice experimental detonations conducted in 1989 
during the MITAVHOI CEAREX ice camp [20], indicating that the appearance of SH 
waves in the PRUDEX data may not be an isolated phenomenon. 
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Chapter 4 

Source Location 

As discussed in Chapter 3, a critical part of the analysis of propagation data 
obtained at the PRUDEX ice camp is the two-dimensional localization of the under-ice 
explosive shots which excite seismo-acoustic propagation in the ice cover. Although the 
depth for these shots is known, only the most general information on their range and 
bearing from array center is available from records of the experiment. 

Prior to commencing the analysis, it was believed that the information available 
from up to sixteen hydrophones could provide the basis for source localization, and that 
geophones at only four locations could contribute little of additional use. This 
predisposition to rely on the hydrophone data proved to be entirely erroneous, not only 
because analysis using the hydrophone data was insufficiently accurate, but also because 
it concealed information crucial to a thorough understanding of propagation in the ice. 

This chapter reviews the data available to the localization effort and the routines 
used in localization, then examines the results for both hydrophone and geophone data. 
A detailed analysis of apparent anomalies observed in the geophone results leads to 
some unexpected conclusions about what these anomalies reveal about propagation in 
the PRUDEX ice cover. Finally, possible locations for the source of the SH waves 
noted in Chapter 3 are reviewed. 

4.1 Localization Data 
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To take full advantage of the unique shock wave/bubble pulse character of both 
the received hydrophone and geophone data, the occurrence of each pulse series is 
identified in each received time series and associated with the appropriate non-dispersive 
wave type, i.e., the longitudinal wave, the SH wave or the waterborne pulse response 
in the geophone data, or the waterborne wave in the hydrophone data. Times of each 
pulse maximum and minimum are interpolated to the nearest 0.1 msec. Thus, for each 
geophone time series the basic localization data consist of times of received maxima and 
minima at each of four geophone locations for the shock wave and the first several 
bubble pulses associated with each of three separate wave types. For each hydrophone 
time series, the localization data consist of the times of received maxima and minima 
for the shock wave and the first several bubble pulses at thirteen hydrophone locations 
(hydrophones 6, 9 and 12 in Figure 3-1 were not connected to the recording system). 

The accuracy of the localization data so developed is effected by the limitations 
of the recording system and the interpolation routine, as well as by the presence, 
principally in the geophone data, of interfering waves. In order to estimate the 
uncertainty in the localization data, the hydrophone maxima are averaged to determine 
the mean time difference between a shot’s shock wave and each of its bubble pulses. 
Since this bubble pulse At should be some constant value regardless of wave or receiver 
type, deviations from the mean value are used to develop the statistics of both 
hydrophone and geophone data so that an appropriate weight can be assigned in the 
localization routine. 

A second key element of localization data is the sound speed in the first few 
meters of sea water directly beneath the ice. Fortunately, this data is readily available 
[16] for times within one hour of the shot times of interest. Figure 4-1 is a typical plot 
of temperature, salinity, and sound speed for the PRUDEX ice camp. 
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Figure 4-1: Temperature/Salinity/Sound Speed profiles at the PRUDEX ice 
camp, 31 March 1987, 0601 U.T. (from McPhee [16]). 

Although not a major factor in the localization, the depth of each shot was 
recorded during testing. Informal discussions with personnel involved in the experiment 
indicate that these depth measurements were obtained by lowering the charges into the 
water on a marked line and are reliable values. 

4.2 Localization Routine 

To provide a flexible tool with which to localize each shot in the specialized ice 
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geometry, a FORTRAN computer routine centered around a singular value 
decomposition (SVD) routine adapted from Press, et. a/.[21], is employed. SVD 
analysis provides a consistent method of combining data points from many different 
array elements, shots, and wave forms, while maintaining the ability to supply 
supporting parameters which might be known separately, such as wave speed and shot 
time, and the ability to specify the uncertainty in all data points and parameters 
individually. A system of equations is assembled in matrix form based on a simple 
linear relation for each pulse arrival at each array element, 

( 4 . 1 ) 

where tj‘ is the shot time for the jth shot, the received time for a specific pulse of 
wave w at array element a from the jth shot, r, is the assumed range from the shot to 
array element a (treated as a known value), and 1/v^ is the inverse speed of wave w. 
For example, a simple matrix system consisting of the shock waves from two different 
wave types (e.g., the longitudinal wave and the SH wave), arriving from two separate 
shots conducted at a single location, and received at a two element array, becomes 
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(4.3) 



If one of the speeds or shot times is known, an additional row is appended to matrix A 
and vector b with this information. The system is then weighted by dividing each row 
of matrix A and vector b by the uncertainty of that row (making the reasonable 
assumption that the uncertainty in each equation is uncorrelated). For a given shot 
location, singular value decomposition solution of (4.2) after weighting yields the best 
fit to the data for shot time and wave speed. The uncertainty of each estimate is also 
available [21]. 

In order to locate a shot in space and time, the above SVD routine is employed 
in a matched field approach to search through two dimensional space for a best fit for 
shot times and wave velocities, as indicated by the lowest value of chi-square (x^), the 
sum of the square of the difference between the modelled array times and the data. The 
best fit corresponds to the most likely location for the shot. To evaluate the uncertainty 
in the best fit location, values of x^ are determined for each point in space searched by 
the localization routine, converted to dB, normalized to zero dB at the lowest value of 
X^ (i.e., converted to curves of ax^ from the best fit), and plotted. The best fit location 
and the estimated uncertainties (assumed to be Gaussian) are then used to compute a 
large number of Monte Carlo simulations of the data [21]. These Monte Carlo 
simulations are supplied to the localization routine and their best fit locations 
determined. Simultaneously plotting these best fit Monte Carlo-derived locations with 
the AX^ contours resulting from the original best fit allows a straightforward assignment 
of confidence limits to specific dB values of ax^. 

4.3 Localization with Hydrophone Data 
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Although use of all available elements in the hydrophone array was intended, 
when localization using hydrophone data was initially attempted, it became immediately 
apparent that positions recorded for the elements located 100 meters and further from 
the array center were very inaccurate. Localization failed to converge meaningfully 
when data from these outer hydrophones were incorporated, and attempts to refine outer 
hydrophone positions with data from different shot locations on successive days failed 
due the small angular separation and long time interval between shot series. Limited 
thus to the inner 9 hydrophones (of which 8 were connected to the recording system) 
the localization routine proved very successful in refining shot bearing, but much less 
so in refining shot range. 

Figure 4-2 is a plot of contours for a best fit shot location determined using 
hydrophone data, along with the best fit locations for 80 Monte Carlo simulations of that 
data. The 1 dB ellipse, corresponding to about a 90% confidence limit, has a major axis 
more than 900 meters long. With an estimated range from shot to array center of much 
less than this distance, ^Yave velocities determined using travel times from the best fit 
shot location cannot be specified even to within ±50%. As will be seen in the next 
chapter, inversion of ice propagation data is heavily dependent on measured wave 
velocities. Clearly, hydrophone-based shot locations do not provide the accuracy which 
reliable inversion requires. 

4.4 Shot Location using Geophone Data 

After having failed to determine a reliable location and time for any experimental 
shot using the hydrophone data, localization was begun using geophone data. With the 
geophone data it was hoped that the information carried in the multiple waves could 
offset the limitations of an array comprised of only four more closely spaced element 
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Figure 4-2: Contours for the best fit location for the F series of shots 

calculated using hydrophone array data, plotted with the best fit locations(+) for 
Monte Carlo simulations of that data. 

positions. The results eventually much more than justified the expectations, but not until 
a thorough analysis and explanation of some apparent anomalies was completed. 

4.4.1 Variations in Ice Thickness at the Receiving Array 

The first anomaly is relatively easy to understand and eliminate. When only the 
arrivals of the response to the waterborne pulse are processed in the localization routine, 
the shot bearings correspond well with those produced by the hydrophone data, as is 
expected, but wave velocities are consistently reported as 10-15 m/s higher than the 
known value of about 1435 m/s. Since the relatively low uncertainty of the pulse 
measurements supports more accurate velocity determinations, some variation in the 
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travel path is postulated. Detailed numerical studies show that this problem is very 
likely due to variations in the ice thickness below each geophone - changes in pulse 
travel time of as little as 0.3 ms, as would be caused by a difference of 1 meter in ice 
thickness beneath different geophones, account for the velocity difference without 
significantly altering the bearing reported by the localization routine. To estimate the 
thickness at each geophone, ice thicknesses at the geophones are added as variables to 
the SVD system in the localization routine, so that (4.1) becomes 






(l 









(4.4) 



where a is the (known) compressional wave velocity and 2h, the unknown thickness at 
geophone a. To anchor the estimates, since the travel times are only weakly dependent 
upon ice thickness, thickness variables are also added as parameters and assigned the 
same average value and expected error. With this modification, the localization routine 
is applied to the existing best fit location and values for the geophone ice thicknesses 
determined. These values are then used to correct the received times, and the basic 
localization routine is used to determine a new best fit location. This location can be 
used recursively to determine new estimates of ice thickness and improve the best 
position; however, the process generally converges very rapidly. Final values for ice 
thickness at each geophone are shown in Table 4-1. Because the system is very 
sensitive to differences in thickness, but not to absolute thickness, only the differences 
from some average value are shown. 



4.4.2 Evidence for Refracted Waves 
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Geophone 


Deviation from 


Nr 


average thickness 


1 


+0.08 m 


2 


-0.18 m 


3 


-0.83 m 


4 


+0.97 m 



Table 4-1: Deviations from the average ice thickness determined at each 
geophone in the PRUDEX airay. 

The basic version of the shot location routine utilizes only the two most clearly 
defined and essentially non-dispersive wave pulse arrivals at the geophone array, the 
longitudinal wave and the response to the impact of the waterborne pulse on the ice 
under the geophone. Inherent in the basic location routine is the assumption that the 
water and the ice plate are homogenous media, such that non-dispersive waves travel 
in straight lines at constant speed. As seen in Figure 4-1, this assumption is valid for 
the waterborne anivals (except as noted above). For longitudinal waves traveling in the 
ice, however, the difficulties introduced by variations in ice thicknesses at the array 
foreshadow the final conclusion that the PRUDEX ice cover cannot be treated as a 
single homogeneous plate. When both waterborne and longitudinal waves are processed 
in the location routine, the resulting values are much higher than the apparent 
measurement uncertainty supports; however, most of this error lies in very poorly 
predicted arrival times for the longitudinal wave. Detailed review of solutions for all 
eight shots from four different locations discloses that in all cases the longitudinal wave 
discrepancy arises because the longitudinal waves are arriving from an apparent bearing 
about 8 to 10 degrees to the right of the waterborne pulses. A further review for any 
single shot location of the arrivals from all three wave types propagating in the ice 
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shows that each wave type is coming from a different apparent bearing than the 
waterborne pulse. Figure 4-3 is a plot of the apparent direction of wave arrival 
computed for each wave type from a single shot location. Assuming that all waves 
originate from the same point. Figure 4-3 strongly suggests a family of waves refracted 
in the horizontal plane. 

Experimenters returning from the PRUDEX ice camp described a relatively 
straight ridge in the ice cover which marked the transition from the thin runway ice 
under which the experimental shots were made, to the thicker ice on which the receiving 
array was situated. This description, combined with the indications of refraction noted 
above, prompted an investigation into the possible existence of the horizontal refraction 
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Figure 4-3: Plot of PRUDEX array layout showing apparent axis of arrival of 
(1) both the hydrophone and geophone water waves, (2) the SH wave, (3) the 
longitudinal wave, and (4) the flexural wave. 
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of waves propagating in the ice along a line created by a vertical plane separating the 
ice cover into two homogeneous half-plates of different thicknesses. 



To investigate horizontal refraction, the location routine is further modified to 
determine, for a given index of refraction and ridge line orientation, the path which a 
wave travels from source to receiver. For refracted arrivals (4.1) becomes 
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(4.5) 



where r,. is the distance traveled in the first half-plate from the source to the ridge line 
on a path to array element a, rj, is the distance traveled in the second half-plate from 
the ridge line to array element a, and and the corresponding wave velocities 

in the two half-plates. Finally, adding an additional equation to the system. 



J ^L=o 



(4.6) 



where n is the index of refraction, establishes the appropriate ratio between the two 
velocities. In order to locate the most likely ridge line, a comprehensive search is 
conducted by computing for each shot the best fit shot location and corresponding value 
of for a broad range of possible ridge line orientations and index of refractions. The 
most likely orientation is determined by combining the values for all shots and 
choosing the orientation that produces the lowest for all shots together. Figure 4-4 
shows the contours for various values of the slope and y-intercept of the ridge line 
referenced to the array as an x/y plane with the origin at the center, and the positive x 
axis along the line containing geophone 4 and hydrophones 4, 8 and 12. Figure 4-5 
shows the orientation of the best ridge line, described by the line 

y = -1.95x + 223 meters (4*7) 
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Figure 4-4: contours (in dB) for the best fit ridge line orientation, 

described by the line y=mx+b on an x/y plane centered on the horizontal plane 
of the array with geophone #4 on the x-axis. 
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Figure 4-5: Plan showing on an x-y plane the PRUDEX geophone(*) and 
hydrophone(+) array, the best location for the F series of shots(x), the best fit 
ridge line(— ), and the longitudinal wave paths from shot to each geophone(-). 
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Figure 4-6: Aerial photograph of the PRUDEX ice camp and array, showing 
locations of identifiable hydrophones (geophones and some hydrophones are not 
visible), the array axes and the ice ridge line. 
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Figure 4-7: ax^ contours (in dB) for the best fit location for the F series of 
shots calculated using geophone data, plotted with the best fit locations for 80 
Monte Carlo simulations(+) of that data. 

relative to the array and the F series shot location. For any given shot series, the 
value for the best fit refracted path is fully 6 dB better (lower) than the value for the 
best fit to the same data on an unrefracted path. Figure 4-6 is the only available aerial 
photograph of the PRUDEX camp showing both the ridge line and the array layout. 
The best fit ridge line appears to correspond remarkably well with the ridge line in the 
photograph. 

The search for the best ridge line necessarily includes as part of its operation the 
best fit location for all shots used in the search. Figure 4-7 is a plot similar to 
Figxire 4-2 showing the a%^ contours for the best fit location for the F series shots, as 
well as the results of best fit searches of 80 Monte Carlo simulations of the received 
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data. The major axis of the 90% confidence limit ellipse (corresponding to the 1.6 dB 
contour) is only 60 meters long, implying an accuracy of about ±5% in wave velocities 
determined using this position. In fact, the slower velocities associated with the flexural 
wave will tend to be much more accurate than ±5%, since the localization routine 
determines a corresponding shot time for which a wave traveling with the speed of 
sound in water (1435 m/s for these shots) will be measured exactly, regardless of the 
error in range. As a wave speed increases or decreases from this value, the speed 
measurement error increases accordingly. A wave traveling at 1000 m/s, for instance, 
will be measured at a range of 6(X) ±30 m to within ±1.6%, while a slower wave 
traveling at 5(X) m/s will be measured to within ±3%. 

Figure 4-8 provides a direct comparison of 90% confidence limit ellipses for the 
hydrophone and geophone based best fit locations for a shot series. While a system to 
monitor hydrophone positioning, a more sophisticated processing system, and a larger 
array can certainly improve hydrophone performance, it is remarkable that with all of 
these improvements, the performance of a hydrophone-based system will be unlikely to 
surpass that of a simple system of four 3-axis geophones. 

4.5 Locating the Source of the SH Wave 

The SH wave arrivals are investigated using the shot location routine in an 
attempt to determine their source and time of origin. Analyzing the SH wave arrivals 
separately, the range resolution is, as could be expected from the hydrophone results, 
very poor. Assuming the point of origin of the SH wave to be any given position 
between the shot and the ridge line, however, the best fit time of SH wave origin 
calculated by the location routine is roughly consistent with travel from the shot location 
to that point of origin at about the same speed as the SH wave’s travel from the point 
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Figure 4-8: Plots of the 90% confidence limit ellipses for the F series shot 
location derived from geophone and hydrophone data. 



of origin to the array; i.e., the SH wave arrivals are consistent with generation by the 
shot itself, or with generation by the interaction of some wave traveling at about the 
same speed as the SH wave at some location near a line between the shot and the array. 

Figure 4-9 is a plot of a%^ contours for the SH wave arrivals determined by 
restricting the time of origin to the shot time determined with the procedures of Section 
4.4. The intersection of the 1 dB contours for the two locations demonstrates 
qualitatively the compatibility of the two solutions, and supports the argument that the 
SH wave is excited directly by the shot itself. 

A second possible source for the SH wave is out-of-plane scattering during the 
interaction of the longitudinal wave with the ridge line. Figure 4-10 is a plot of 
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contours for the SH wave arrivals calculated by constraining the time of origin to the 
average arrival time of the longitudinal waves at the ridge line en route from the shot 
to the array. The large offset of the center of the ellipse from the ridge line indicates 
that the SH wave arrivals are not consistent with their generation during the interaction 
of the longitudinal wave at the ridge line. Similarly, the flexural wave interacting at the 
ridge line is another possible mechanism for out-of-plane scattering and production of 
the SH wave; however, that possibility can be immediately dismissed as the SH waves 
arrive at the receiving array before the flexural waves (at the peak frequency of 20 Hz) 
arrive at the ridge line. 

A fourth possible source is an interaction of the waterborne wave with some 




X range (m) 

Figure 4-9: contours (in dB) generated by the shot location routine for the 

SH wave point of origin(-) assuming the time of origin is fixed at shot time, 
with contours for the best fit shot location(-*) of section 4.4. 
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Figure 4-10: contours (in dB) generated by the shot location routine for the 

point of origin of the SH wave with time of origin fixed at the average time of 
longitudinal waves’ (•••) arrival at the ridge line (— ). 

discontinuity on the underside of the ice. As the wave speed in water is comparable 
with that of the SH wave, this possibility cannot be ruled out immediately; however, 
Figure 4-11, the contours for the SH wave point of origin assuming that the SH 
wave time of origin coincides with the arrival of the waterborne wave at the ridge line, 
shows at least that the water wave/ridge line interaction is probably not the source, again 
because the location contours for the time of that interaction are offset relatively far 
fi'om the ridge line. 

Based on available evidence, the source of the SH wave cannot be positively 
identified. Origin at or near the shot location and time is consistent with the data. 
Interactions of the flexural wave, the longitudinal wave or the waterborne wave with the 
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Figure 4-11: a%^ contours (in dB) generated by the location routine for the 
point of origin of the SH wave with the time origin fixed at the time of the 
waterborne waves ’(•••) arrival at the ridge line(— ). 

ridge line do not appear to be likely sources for the SH wave. Since the major known 
discontinuities in the PRUDEX ice canopy are associated with this ridge line, it is most 
likely that the SH wave is either excited directly by the under-ice detonation in an as 
yet unexplained manner, or it is generated by out-of-plane scattering in the immediate 
vicinity of the detonation, perhaps during interactions with some unknown feature on 
the underside of the ice canopy, in which case there is insufficient information in the 
PRUDEX experiment to determine which wave is the source. 
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Chapter 5 

Inversion of Propagation Data 



Solving the wave propagation problem in the manner of Section 2.1, computing 
particle motion based on the equations of motion, the characteristics of the material, a 
given geometry and a known excitation, is the forward problem in seismo-acoustics. 
Although solution of the forward problem is seldom easy, the approach is at least 
straightforward and the correct solution should be unique [7]. Taking the measured 
particle motions in response to a known source and processing that data "backwards" 
through the appropriate equations to obtain the unknown elastic and geometric 
parameters is the inverse problem. Inversion of seismo-acoustic data often becomes 
more complex than the comparable forward problem because it assumes solution of the 
forward problem as a starting point and must deal with non-unique solutions. A given 
set of elastic/geometric parameters will produce only one response to a given excitation; 
however, it is very possible that different sets of those parameters will produce the same 
measured response to that excitation. The likelihood of non-unique solutions to the 
inversion process necessarily increases as the number of unknown parameters increases. 

5.1 Inversion Parameters 

The principal parameters desired from the seismo-acoustic inversion can be seen 
by inspection of the equations of motion in a linearly elastic solid (2.1). The motion 
of ice particles, hence the propagation of waves in the ice, is dependent upon the ice 
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density p, and the ice Lame constants p and X. Equivalently, using relations (2.6) and 
(2.7), the Lame constants may be expressed as the compressional and shear velocities 
a and p. 

Although not expressed directly in the equations of motion, internal friction in 
a propagating medium dissipates the energy of waves propagating in that medium. In 
most cases this attenuation must be known a priori or added to the list of parameters 
to be determined in the inversion. In this work attenuation is described by the two 
parameters y„ and Yp, the attenuation of compressional and shear waves, respectively, 
as described in Section 2.2.2. 

Ideally, there are no unknown geometric factors to complicate an inversion. In 
this chapter shot location relative to the receiving array, as determined in Chapter 4, is 
considered a known value. The sparse information available about ice thickness in the 
vicinity of the PRUDEX ice camp necessitates treating ice thickness as an unknown and 
including it in the inversion. 

5.2 Previous Measurements 

Although the mechanical properties of sea ice have been extensively studied 
[22], very little work has been done to determine the low frequency elastic properties 
of the arctic ice cover in situ. Until recently, actual measurements have been limited 
to some early wave speed measurements in freshwater lake ice [5] [6] [23], and pack 
ice [24] [25], high frequency attenuation measurements in glacial ice [26] [27] 
and sea ice [28], wave speed profiling in both lake and sea ice [29], and data 
obtained fi*om small scale laboratory experiments. As a result, determinations of the low 
frequency properties of arctic sea ice were largely inferred from other ice environments 
or extrapolated from high frequency laboratory and in situ data. As an excellent 
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example of this approach, McCammon and McDaniel [22] have employed a 
comprehensive summary of available laboratory and field measurements to determine 
values for the attenuation in the arctic sea ice, for use in studies of the acoustic 
reflectivity of the ice cover. Based on this summary, they have estimated that 
compressional wave attenuation can be approximated by 

Y. = eio-’o ^ (5.1) 

X 

and (assuming Poisson’s ratio to be constant at 0.33) shear wave attenuation by 

Yp = 3.610-^P ^ . (5.2) 

Results obtained by such a combination of extrapolation and inference may 
accurately represent the characteristics of seismo-acoustic propagation at high 
ft’equencies, but it is questionable whether these values may be translated to reflect low 
frequency behavior as well. At low frequencies and long wavelengths, macroscopic 
discontinuities, such as cracks and ridges in an arctic ice plate, may reduce propagation 
speed and increase attenuation in the medium. 

Recently, several investigators have obtained values for the elastic parameters of 
arctic ice at low frequencies. In 1986 Stein [14][30], as an adjunct to other studies, 
estimated values for shear and compressional velocities and attenuations from earlier 
work at two arctic sites. In 1989 Brooke and Ozard [31] completed a detailed study 
of the elastic properties of sea ice based on measurements in the Slidre Fjord of the 
Canadian Archipelago in 1986 and 1987. The results of Stein, and Brooke and Ozard 
are summarized in Table 5-1. 



5.3 The Inversion Procedure 
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Stein 




Brooke and Ozard 
Smooth Ice 


Rough Ice 


data date 


1980-82 




1987 


1986 


1987 


1986 


a m/s 


3500 




NA 


NA 


NA 


NA 


CipITl/S 


NA 




3084 


2960 


2893 


2864 


P m/s 


1800 




1705 


1891 


1660 


1746 


y^dB/k 


0.46 




NA 


NA 


NA 


NA 


7p dBA 


1.57 


20-40HZ 


0.32 


0.45 


2.33 


1.26 






40-80HZ 


1.00 


0.57 


2.55 


0.84 






80-120HZ 


0.38 


0.49 


1.33 


0.48 



Table 5-1: Summary of recent measurements of the elastic parameters of arctic sea 
ice at low frequency. 



Inversion of the PRUDEX propagation data is conducted initially with the 
assumption that the floating ice canopy between source and receiver can be treated as 
a single homogeneous plate, and that only the longitudinal and flexural waves are 
excited by the explosive charge. Although neither of these assumptions is actually valid, 
the methods of this procedure serve to demonstrate the power of SAFARI modeling, and 
the results can be viewed as a form of "average" behavior. The inversion is then revised 
to reflect the more complete knowledge of the plate’s character obtained in Chapter 4. 



5.3.1 Inversion for an Infinite, Homogeneous Plate 

A straightforward way to simplify a seismo-acoustic inversion is to select a small 
subset of the full set of elastic parameters describing the propagating media, and isolate 
for study a portion of the measured response which is sensitive only to the elements of 
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this subset Preliminary study of the sensitivity of the flexural wave in a floating ice 
plate to variations in the elastic parameters shows that the flexural wave is relatively 
insensitive to variations in both the compressional wave velocity and the density of the 
ice. Because of the insensitivity to the density of ice, a nominal value of .91 gm/cm^ 
[32] is used throughout this work, and no attempt is made to determine ice density 
in any inversion. More importantly, this study shows that the flexural wave is 
dependent only upon the shear velocity, the attenuation values, and the thickness of the 
ice - it can be isolated to invert only for this more limited number of parameters. In 
addition, the dispersion curve of the flexural wave, i.e., the relation of the flexural 
wave’s group velocity to its frequency, is essentially independent of attenuation, and the 
measured dispersion curve can be inverted for only the shear velocity and the ice 
thickness. 

Adopting an approach similar to that used by Jensen and Schmidt [33] to 
determine shear speed and shear attenuation of the sea bed from analogous Scholte wave 
data, the first step in the inversion consists of constructing a dispersion curve for the 
flexural wave. The dispersion curve is developed using the flexural wave responses and 
the positions determined for the eight experimental shots which not only excited a 
vigorous flexural wave, but also excited observable longitudinal waves and responses 
to the waterborne pulse (necessary for localization). In this section the positions used 
are those determined by treating the ice sheet as a homogeneous plate. The dispersion 
curve is built by applying a moving Fourier transform with a Hanning window to the 
time series, with the necessary window size determined by the simple expedient of 
increasing its length until the dispersion curve is stabilized [34]. All thirty-two 
individual curves (eight shots received at four vertical geophones) are then normalized 
to a constant noise value and combined to produce a single best dispersion curve for the 
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assumed homogeneous ice plate. 

Having obtained an experimental dispersion curve, the inversion proceeds by 
selecting a likely model for the unknown parameters, using the SAFARI numerical 
modeling routine to determine a synthetic dispersion curve, comparing the model and 
the experimental dispersion curves, developing a correction to the model, and recursively 
refining the model until the model’s curve converges with the observed one. If the ice 
thickness is known, the above procedure will quickly determine the correct shear 
velocity; unfortunately, the ice thickness at the PRUDEX ice camp is not known. For 
any given thickness and shear velocity, a family of solutions are found which can 
reproduce the same flexural wave velocity at a given frequency simply by adjusting the 
ice thickness up or down and compensating with an appropriate change in shear 
velocity. While some difference does arise between two such similar solutions over the 
frequency range of interest (2-60Hz), this difference is well within the accuracy 
available in comparing dispersion curves. This problem is illustrated in Figure 5-1, 
which shows several families of dispersion curves calculated for two different shear 
velocities and various ice thickness values. 

The uncertainty in shear velocity is largely eliminated by expanding the inversion 
to include the longitudinal wave. As the longitudinal wave is essentially non-dispersive, 
this expansion necessitates a shift to the time domain and a direct comparison of 
experimental and synthetic time series. In addition to providing the compressional wave 
velocity, this expansion has the added benefit of allowing the estimation of the 
compressional and shear wave attenuation as well. The longitudinal wave is largely 
insensitive to ice thickness, casting some doubt on its ability to diminish the uncertainty 
in the shear velocity; however, it is so sensitive to shear velocity that only a very 
limited range of shear velocities can combine with a reasonable compressional wave 
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Figure 5-1: Two sets of dispersion curves for flexural waves in ice at shear 
velocities P=1600m/s(--) and P=1800m/s( -)> and (top to bottom in each set) ice 
thicknesses of 1.25, 1.15 and 1.05m. 

velocity to match the observed longitudinal wave. Figure 5-2 illustrates the dependence 
of the longitudinal wave on shear and compressional velocities, as well as demonstrating 
its essentially non-dispersive nature. Of equal importance in reducing the uncertainty 
in the uniqueness of the inversion, matching the observed flexural wave in the time 
domain is more sensitive to errors in flexural velocity over a wide frequency range than 
matching the calculated and observed dispersion curves. 

A complication introduced by the decision to shift to comparing synthetic and 
observed time series is the necessity to provide an accurate representation of the 
explosive source used to excite the observed waves. A computer routine based on 
equations provided by Wakeley [35] has proven to be very successful in reproducing 
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the acoustic pressure signature measured at 1 meter from a known underwater explosive 
source; although possibly due to the effect of the extreme cold on the explosives, the 
bubble pulse intervals predicted by Wakeley’s routines are consistently longer than 
observed at PRUDEX for the same explosive weight and depth. This discrepancy is 
resolved by reducing either the explosive weight or the proportionality constant used in 
the equations slightly from that provided by explosive tables [36] [37] for a given 
length of primacord or dram weight of explosive, such that the synthetic and observed 
bubble pulse intervals agree. In this way the relative spectrum levels of the real and 
synthetic shots are identical. It is likely that the absolute levels are also equal, but there 
is insufficient information available to verify this assumption. Figure 5-3 shows the 




Figure 5-2: Two sets of group velocity curves for longitudinal waves in ice 
with compressional velocities of 3500m/s(-) and 3400m/s(— ) and shear 
velocities (top to bottom in each set) of 1800, 1700 and 1600 m/s. 
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Figure 5-3: Synthetic time series for the pressure signature at 1 meter for an 
explosive charge simulating shot F3; tO£, sampled at 10 KHz, and bottom. 
prefiltered and decimated to 1000 Hz. 





Figure 5-4: Top, synthetic time series of Figure 5-3 filtered to a 2-90 Hz band, 
and bottom, spectrum of filtered time series. 
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synthetic pressure signature of the explosive charge (1 foot of primacord) used in shot 
3 of series F. Note that after prefiltering and decimation to a sample rate of 1000 Hz 
the first bubble pulse is larger in both amplitude and energy content than the initial 
shock wave, duplicating the relation seen in the experimental pulse trains. One of the 
requirements of the SAFARI pulse calculation routine is that to avoid "ringing" in the 
output time series, the frequency integration routine must be truncated where the source 
pulse has a frequency minimum [3]. To meet this requirement, as well as to limit the 
computation required for the frequency integration, the source pulse is digitally filtered 
[38] [39] to a 2-90 Hz band, as shown in Figure 5-4. Limiting the frequency 
integration to this band has no effect on the inversion. The partial spectrum of a typical 




Figure 5-5: Spectrum of signal received on vertical component of geophone #3 
during experimental shot F3, showing preponderance of energy in the 2-90 Hz 
band. 
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geophone time series. Figure 5-5, shows that most of the information carried in the 
signal resides in that 2-90 Hz band. 

Figure 5-6 shows the measured dispersion contours and Figure 5-7 the dispersion 
contours determined from a best fit synthetic time series, both calculated for the best fit 
shot location determined by treating the PRUDEX ice cover as a single homogeneous 
plate. Both figures also include the exact dispersion curve calculated for the 
homogeneous plate’s best fit inversion parameters. Two figures which demonstrate the 
power of SAFARI pulse modeling are the plots of synthetic and observed geophone time 
series for shot number 3 of series F, Figure 5-8 for the horizontal geophone and 
Figure 5-9 for the vertical geophone. Note in both figures that not only is the flexural 
wave modelled well, but also the longitudinal wave and the response to the waterborne 
acoustic pulse. Other arrivals seen after the longitudinal wave, but before the flexural 
wave are probably due to the inhomogeneity/anisotropy of the real ice and are not 
reflected in SAFARI modeling. Also of interest, the apparent irregular behavior on the 
tail of the synthetic flexural wave is introduced when the air is modeled with a realistic 
sound speed and density rather than treated as a vacuum; however, no set of air 
parameters modeled the real response well, and the air is treated as a vacuum for the 
remainder of this chapter. 

5.3.2 Inversion of Two Abutting Infinite Half-Plates 

As discussed in the previous chapter. The PRUDEX ice plate is much more 
accurately described as two half-plates of different thicknesses. Although a version of 
SAFARI able to handle some range-dependence, including inclusions in an ice plate of 
a different thickness than the rest of the plate, is under development during the summer 
of 1990 by Gerstoft and Schmidt [40], it is not available as of this writing. In order 
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frequency (Hz) 

Figure 5-6: Observed contours of spectrum level (in dB normalized to 0 dB 
maximum) obtained by combining data from 8 shots at the PRUDEX ice camp, 
with the dispersion curve (-•) calculated for P=1700m/s and 2h=1.31. 
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Figure 5-7: Synthetic contours of spectrum level (in dB, normalized to 0 dB 
maximum) derived from SAFARI time series calculated for a=3400 m/s, 
P=1700 m/s, 2h=1.31m, with corresponding exact dispersion curve (— ). 
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Figure 5-8: Observed response(-) on radial geophone #3 for PRUDEX shot F3; 
SAFARI synthetic radial response(--) for a=3400m/s, P=1700m/s, 2h=1.31m, 
Y„=1.0dBA, Yp=2.99dBA. 
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Figure 5-9: Observed response(-) on vertical geophone #3 for PRUDEX shot F3; 
SAFARI synthetic vertical responsef--) for a=3400m/s, p=1700m/s, 2h=1.31in, 
Y„=ldBA, Yp=2.99dBA. 
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to provide a rough inversion simulating the environment at the PRUDEX ice camp, a 
modification of the stationary phase method [4] is used in conjunction with range 
independent SAFARI solutions to approximate the wave form received in the second 
half-plate after detonation in the first half-plate. 

As seen in Chapter 2, the flexural wave in an ice plate is a strongly dispersive 
wave of a single mode. Following the development in Aki and Richards [4], a wave 
packet composed of a single mode may be expressed as 

. <5.3) 



where |F(co)| is the spectral density and (t>(CD) the initial phase. In stationary phase 
analysis, the integration path is along the real to axis, and for large values of x and t the 
integrand cot +kj^x oscillates rapidly, with each oscillation tending to cancel the next in 
the integral. Only at or near a saddle point, given by 

—(-(jjt+kx)=0 (5.4) 

du> 

will the phase vary slowly enough to provide a significant contribution to the integral. 
(5.4) can be simplified to 



£ 

t 



du> 




(5.5) 



where u is the group velocity. Solving (5.5) yields <Os(x,t), the frequency expected to 
dominate at distance x and time t. Expanding the phase -cot + in a Taylor series 
about the point C0=C0s, and neglecting higher order terms, yields 



-ci>t+k^ “ -o)^r+kj(o)^+ 



xfk. 

2 da 






(5.6) 
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Substituting (5.6) into (5.3), and simplifying, results in 






2n 



^ 2ti 



2 



dX 









(5.7) 



where ± corresponds to d^kydO)^ <0 or >0, respectively. 

In order to apply this method to the problem of propagation in two half-plates, 
the stationary phase approach must be expanded to account for two propagating media. 
If the integrand in (5.3) is modified to reflect propagation for distances x, and Xj in 
media with horizontal wave numbers and k^ 2 > then (5.4) becomes 

= 0 (5.8) 

or 






= t 



for which a solution is 



(5.9) 






du> 



d(a 



(5.10) 



Expanding the two relations in (5.9) in a Taylor series and adding the results yields 



-0(t,+f2)+Vl-^^x2^2 “ 

- « 5(^1 + ^2) + >1 A + 



Xj X2 dX 



x2 



2 do>^ 2 



(w-o^^ 



(5.11) 



Substitution in (5.3) gives 
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(5.12) 



(5.12) does not involve attenuation, so to provide an ability to compare synthetic curves 
with experimental ones more directly, an additional attenuation factor is used: 



/a = ^ 



-(6a^aX*6pt^x*Ax+B) 

f 



(5.13) 



where 5^p=Ya,p*(loge/40jt), x is the distance between the point at which the pulse 
spectral density is determined and the point for which the pulse is being calculated, and 
A and B are constants determined empirically by comparison with known results. 

To demonstrate the potential accuracy of the stationary phase approach applied 
to the flexural wave in ice, SAFARI is used to generate synthetic time series and phase 
and group velocity curves for the response to an explosive shot for a given set of elastic 
parameters in a single homogeneous infinite floating ice plate at a ranges of 242 m and 
569.5 m. The spectral density of the short range shot and the phase and group velocity 
curves are supplied to a computer routine which uses the stationary phase approximation 
(5.7) to generate the curve at the longer range. Figure 5-10 shows the short range time 
series, and Figure 5-11 shows the longer range time series as generated by SAFARI and 
as calculated from the shorter range pulse using stationary phase and the empirical 
amplitude attenuation of (5.13). Although the stationary phase result is not perfect, it 
is nonetheless good enough to allow at least an estimation of best fit curves. 

To apply modified stationary phase to the problem of two abutting half-plates, 
a computer routine is employed which takes the SAFARI-generated phase and group 
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Figure 5-10: Synthetic time series for the flexural wave in a floating ice plate, 
calculated by SAFARI for a=3500m/s, |3=1750m/s, 2h=2.4m, Y„=1.0dBA, 
Yp=2.99dBA at a range of 242m. 
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Figure 5-11: SAFARI synthetic time series(-) for parameters of Figure 5-10 at 
a range of 569.5m, and time series(-*) generated by applying the method of 
stationary phase to Figure 5-10. 
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velocities for the flexural wave in the two half-plates along with the orientation of the 
ridge line separating the two halves and determines as a function of frequency the index 
of refraction and the path from the shot location to a given geophone. The routine uses 
the calculated path, the supplied spectral density of the shot determined from a time 
series generated for the travel distance in the first plate, and the phase and group 
velocity curves to construct using (5.12) the time series as received at the geophone in 
the second plate. This synthetic time series for the flexural wave can now be compared 
with the observed time series to adjust model parameters and proceed with the inversion. 

Reliance on the flexural wave for inversion resurrects the problem of shear 
velocity/ice thickness ambiguity discussed earlier. To establish a reasonable anchor on 
shear velocity, the location routine is used to investigate the best fit SH wave speed in 
the two plates. Assuming that the SH wave is excited directly by the shot and there is 
no anisotropy, this value can be used directly as the shear velocity in the respective 
plates. The best fit consistent with the longitudinal wave’s index of refraction of 1.12 
is Pi=1590 m/s, and p2=1750 m/s, for an index of refraction of 1.10. 

The modified stationary phase procedure introduced above adds an additional 
ambiguity to the inversion problem at any given geophone. A best fit curve can be 
generated by any one of a family of solutions whose flexural wave velocity curves and 
resulting index of refraction combine to produce the same radial velocity from source 
to receiver. If the second plate is truly a homogeneous half-plate of constant thickness, 
then this ambiguity can be resolved with the straightforward but laborious procedure of 
locating the best fit simultaneously at all four geophone locations; unfortunately, 
apparent variations in the second plate’s thickness have prevented finding any set of 
parameters for the second plate which provide a good fit at all, or even most geophones. 
This problem is not surprising in light of the significant variation in ice thickness 
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calculated at each geophone in Chapter 4, and displayed in Table 4-1. To determine the 
correct solution, the location routine is used to determine the best index of refraction for 
the flexural wave, given the shot location and the orientation of the ridge line. The best 
fit solutions with indexes of refraction centered about this value are chosen, assuming 
a constant thickness in the first plate, but an average thickness varying with geophone 
path in the second. 

Figure 5-12 through Figure 5-15 show the best fit flexural wave time series at 
each of the four vertical geophones, along with the experimentally observed time series. 
Since no geophone calibration data is available, a best fit geophone calibration factor 
of 10'^ m/s/volt is applied to all four geophone outputs to allow comparison of 
calculated and observed time series. Solutions at the four geophones are summarized 
in Table 5-2. Figure 4-5 shows that the ray path from the shot location to geophone #1 





Geophone Nr 


Parameter 


1 


2 


3 


4 


«i (m/s) 


3000 


Pi (nVs) 


1590 


Y„,Yp (dbA) 


1.0, 2.66 


2h, (m) 


1.18 


«2 (ni/s) 


3500 


p2 (m/s) 


1750 


Y„,Yp (dbA) 


1.0, 2.99 


2 h 2 (m) 


2.40 


2.37 


2.15 


2.20 



Table 5-2: Best compressional/shear velocities and attenuations and plate 
thicknesses determined by treating the PRUDEX ice cover as two abutting half- 
plates, with the shot conducted under plate 1, and the receiving array on plate 2. 
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Figure 5-12: Observed flexural wave time series(-) for PRUDEX shot F3 at 
vertical geophone #4, and synthetic time series(-*) for shot F3 at geophone #4 
developed using the parameters of Table 5-2. 
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Figure 5-13: Observed flexural wave time series(-) for PRUDEX shot F3 at 
vertical geophone #3, and synthetic time series(--) for shot F3 at geophone #3 
developed using the parameters of Table 5-2. 
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Figure 5-14: Observed flexural wave time series(-) for PRUDEX shot F3 at 
vertical geophone #1, and synthetic time series(-0 for shot F3 at geophone #1 
developed using the parameters of Table 5-2. 
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Figure 5-15: Observed flexural wave time series(-) for PRUDEX shot F3 at 
vertical geophone #2, and synthetic time series(--) for shot F3 at geophone #2 
developed using the parameters of Table 5-2. 
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nearly coincides with the path to geophone #2, and that the path to geophone #3 follows 
exactly the path to geophone #4. Given that Table 5-2 shows that the average ice 
thickness seen by geophones #2 and #3 is less than that seen by geophones #1 and #4, 
the observation that the best fit thickness at geophones #2 and #3 is less than that for 
geophones #1 and #4 is a confirmation of at least the general validity of theses results. 
Even so, uncertainties inherent in the determination of the shot location and the 
orientation of the ridge line, as well as in the application of the modified stationary 
phase procedure and the inversion process itself, all combine to render the values of 
Table 5-2 as no more than estimates of elastic parameters at the PRUDEX ice camp. 
While Table 5-2 should provide a fair representation of that environment, due to the 
complex interaction of factors involved in their derivation, it is not possible to assign 
definite uncertainties to these parameters. 
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Chapter 6 
Conclusion 



The final chapter summarizes the significant results of Chapters 3, 4 and 5, and 
makes some recommendations for future work. 



6.1 Summary 

Work with the propagation data generated at the PRUDEX ice camp has yielded 
a number of significant findings which contribute directly to the body of knowledge of 
seismo-acoustic propagation in the Arctic Ocean. This work has also highlighted the 
importance of certain tools in the furtherance of that knowledge. 

6.1.1 Elastic Parameters of the Arctic Ice 

The values of bulk compressional and shear wave speeds obtained for the thicker 
multi-year ice at the PRUDEX ice camp, 3500 m/s and 1750 m/s, respectively, compare 
very well with similar values obtained by earlier investigators. It is also interesting to 
note that the shear speed measured in the annual ice, 1590 m/s, is considerably lower 
than in the thicker ice, although not as low as some investigators have predicted [41]. 

The work with the PRUDEX data vigorously supports the assertion that useful 
values for the low frequency elastic parameters of arctic sea ice cannot be obtained from 
laboratory measurements or extrapolations from related data. Attenuation values of 
about 1 dB/>. for the compressional wave, and 3 dBA for the shear wave, as estimated 
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in this work, are slightly higher than but consistent with values reported by the two 
previous in situ studies [31][14]; however, these values are more than four time greater 
than the best available numbers estimated using laboratory and other data [22]. 

6.1.2 Propagation Mechanisms 

Study of the PRUDEX data has revealed the presence of strong horizontally 
polarized transverse (SH) waves propagating in the sea ice canopy as a result of small 
underwater explosive detonations. Since the theory of plate wave propagation has no 
mechanism for coupling SH waves in a plate with acoustic waves in an adjacent liquid, 
these waves are entirely unexpected and as yet unexplained. The PRUDEX data sets 
generally support the contention that these waves originate in the ice canopy at or very 
near the time and horizontal location of the detonation. The data sets do not support SH 
wave generation by out-of-plane scattering during the interaction of either longitudinal, 
flexural or waterborne waves with the ridge line identified in the ice sheet. 

This study also has included the first identification of the horizontal refraction 
of a family of wave types propagating in a sheet of arctic ice. Each of the wave types 
appears to obey simple Snell’s law refraction at the linear abutment between the two 
half-plates which comprise the ice canopy, refracting at angles appropriate to the 
different wave speeds in the two half-plates. 

6.1.3 Analysis Tools 

One of the most useful lessons highlighted during analysis of the PRUDEX data 
was the striking superiority of a simple system of four 3-axis geophones over a system 
of nine hydrophones in a larger array. Not only was the geophone array dramatically 
superior in localizing the underwater detonations which excited elastic waves in the ice. 
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it also allowed the detection and study of wave types and propagation phenomenon not 
visible in the hydrophone data. The ability of the geophone array to isolate several 
different wave types traveling at different speeds placed a much stronger bound on 
possible ranges and bearings from the array center to the shot location than did the 
hydrophone array’s reception of the single waterborne wave, despite the fact that the 
hydrophone array was larger and the pulse arrivals could be measured more accurately 
by the processing system. The detection of the SH waves by the geophone array, as 
well as the strong and clear reception of both longitudinal waves and flexural waves, 
could not have been accomplished from hydrophone data. Indeed, tomographic studies 
of acoustic propagation under the arctic ice have shown that the characteristics of the 
ice cover appear primarily as second order effects (e.g., beam displacement of the 
reflected waterborne pulse) in the hydrophone data [42] [43]. Clearly, a geophone 
array is a superior tool for use in studies of seismo-acoustic propagation in a localized 
section of sea ice. 

Inversion of the PRUDEX data has also served to reemphasize the value of 
SAFARI numeric modeling in seismo-acoustic propagation problems. In a homogeneous 
plate the results of Chapter 5 indicate that SAFARI is capable of fully and accurately 
reproducing all of the elements of the real seismo-acoustic signature: the longitudinal 
wave, the flexural wave, and even the response to the waterborne acoustic wave as it 
passes beneath the geophones. 

Finally, a potentially useful tool to extend two-dimensional SAFARI to a range- 
dependent environment, the modified stationary phase approximation of the flexural 
wave (or any highly dispersive wave), has been demonstrated. This method serves as 
an effective if somewhat limited interim fix for the solution to propagation in two 
adjacent plates until such time as development of the next generation of SAFARI-like 
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algorithms is complete. 



6.2 Future Work 

The great variation in the sea ice elastic parameters determined in the vicinity 
of the PRUDEX ice camp, as well as the general and temporal variability reported 
recently by Brooke and Ozard [31] and earlier by Hunkins [24], all indicate that more 
work in this area is appropriate in order to establish a set of parameters which 
characterize the arctic environment accurately over a given area and season. This thesis 
has shown that in a well-surveyed homogeneous environment {i.e., with all geometric 
uncertainties eliminated), basic SAFARI modeling of under-ice detonations is readily 
capable of yielding very accurate determinations of bulk shear and compressional 
velocities and attenuations from geophone data obtained in situ. 

If obtaining geophone measurements in sufficient number and at enough 
locations to accurately characterize the arctic environment proves to be impractical, 
other approaches, such as ocean acoustic tomography, may also be capable of obtaining 
average values of the elastic parameters over large areas. High frequency cross-hole 
tomography conducted directly in the ice [29] can shed additional light on the variability 
of the anelastic properties of sea ice, although extending such results to the low 
frequencies of interest in this work will remain a problem. These approaches certainly 
warrant further investigation. 

Much more work remains to be done to determine the mechanism which couples 
acoustic waves in the water with SH waves in the plate. An important element in this 
work should be under-ice explosive shots made with geophone detectors installed not 
only at a central array, but spread in range along the propagation path. Additionally, 
under-ice surveys should be conducted to identify discontinuities in the ice canopy not 
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visible on its upper surface. The information available in this expanded experiment may 
well prove vital to isolating the interaction which generates the SH waves. 

In order to further study seismo-acoustic propagation in arctic ice in a range- 
dependent environment, it will be necessary to bring advanced versions of SAFARI, now 
in development [40], to bear on the problem. In this way experimental measurements 
will not necessarily be limited to strictly homogeneous environments, and data sets taken 
in complex environments, such as that of the PRUDEX ice camp, can be inverted with 
more confidence and reliability than is possible with the limited tools now available. 
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